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Integration of ordinary differential equations with 
automatic selection of elementary intervale (DOI or 
SADOI ) 

Closed - one program parameter -~ entered with link in 
Q« Upon being entered, the routine finds an optimum 
elementary interval, executes one step at that interval 
and returns control to the link address* 

Normal entry is at left of (F7) + 1 if (F7) is the location 
of the first word of this routine <> 

Special "initial entry"; To start an integration off, put 
the initial values of the independent variable and. the 
dependent variables into their data bank locations, put 
the appropriate program parameter into location and. 
enter at right of (F7) + 0* The initial sense of integra- 
tion is always such as to increase the independent variable 
x, 

Special "reversing entry"; In order to effect a reversal 
of the sense of integration, enter at right of (F7) + 1 
instead of at lefto Subsequent normal entries will 
cause x to continue changing in the new sense until 
another reversing entry is made., 
210 plus data bank of lOn + 2 words . 
Locations 0, 1, 2j D + 7n to D + lOn - 1 inclusive . 
Locations 3 and k must contain the following parameters 
during input of this routine; 

Significance 
D specifies the location of the data banko See description 
of data bank. 

n is the order of the system of differential equations, 
or the number of dependent variables obeying first order 
equations. Note that in order to revise D and/or n the 
F7 program must be read in over again 8 
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PROGRAM PARAMETER: 



The program parameter must "be supplied only at initial 
entry , when it must be put into location "before enter- 
ing <, It is a compound word of the structure: 



(A - m + l) 2~ 7 +oC 2~ 19 + (e + l) 2* 



■27 



+ X 2 
o 



■39 



Here the integer jJ^ specifies the 

of elementary interval to be permitted! m is an integer 



maximum size h" = 2 
o 



/->■ 



scaling parameter for the derivatives! oC is the location 
of the first word of the auxiliary subroutine! and e is 

an integer specifying the accuracy asked for c The precise 

th 
meaning of e is that the e binary digit in every y. is 

to be correct after integrating from x to x + 1.. 

00 o : o 

The allowable ranges of the integer parameters are: 

< m < 39 
1 < j2 < 39 and m < Jl, 



JL " m < e < 'C =m+39 and e < 62 



o — — "o 
If the integer parameters violate these restrictions 3 \ then 
either the FF stop or a zero-shift stop will occur 
Note that M- 3 m., ^ and e may be revised without reading 
the F7 program in over again,, by making a new initial 
entry * 

DURATION: Normally about (21n milliseconds + 2T^- ) per elementary 

step,, where T> is the time taken by the auxiliary sub- 
routine to calculate derivatives « About 6n milliseconds 
additional is required when the interval must be revised* 
The initial entry., but only the initial entry, requires 2k 
or more times the above time*. 

FF STOP DURING OPERATION: 



Order with 
Sexad. Address 

FF 0^-J 



From 
Location: 

(F7)+0 left 



-39 



Significance: 
The tests indicate that an interval < 2 
would have to be used* Either faulty program 
(e.g.,, overflow) or too large an e for single 
precision operation exists. 
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FORM OF EQUATIONS; In order to use this routine the system of ordinary 

differential equations to "be integrated must be east 
into the form 

£- = f . (x, y v y 2 , - y n )j i = 1, 2 - n 

with x and y. real numbers in machine range and f. real 
numbers such that 2 f . are in machine range* Thus equations 
of higher than first order must be reformulated, by intro- 
ducing auxiliary dependent variables obeying first order 
equations; and the variables of the original scientific 
problem will in general have to be suitably scaled to 
make the problem conform to the above . 
ELEMENTARY INTERVALS WHICH CAN OCCUR: 

The computer's choice of elementary intervals is limited 
to positive and negative inverse powers of 2 (inverse 
jc, to inverse 39 power inclusive) because any 
other choices would much reduce the computing speeds Thus 
if a solution table is to be output with some constant 
interval, the independent variable x must be so scaled, 
that the desired, interval is represented by a power of 
2 in the computer. 
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DATA BAM; 



The data "bank is arranged as follows 



Location 

D-2 
D-l 



Contents 
x 



x 



o 



Significance 

initial value of independent variable 
current value of independent variable 



D y ) 

D+l ,j£ ) 

= „„«=„ = _„-_ = „ = = „-. ) current rounded values of dependent 

) variables 
D+n-1 y n ) 


B+n 2"^ ) 

„ = „-„ = -_.=.____-__ ) current values of scaled derivatives 

D+2n-l 2" m f ) 

n 


D+2n a ^ 2~ m (h/2i) d^/dx 2 ) 

_.= -.-________ = __- ) current scaled. 2nd. derivatives 

D+3n-i a — 2 (h/2i) d yi./dx ) 


D+3n b x £ 2" m (h 2 /30 d^/dx 5 ) 
....__.._,_.,._..___.-.._ ) current scaled 3rd derivatives 

D+lta-1 b = 2~ m (h 2 /3l) d 5 y /dx 5 ) 
n - v ' ' n ; / 


D+ifn c^ 2 =m (h 5 AO dVj/cbr ) 

_ = = .._.._.._...,___ _, _ , ) current scaled. 4th derivatives 
D+5n-l c = 2" D1 (rAi) dy /dx ) 


D+5n d 1 5 2" m (h 751) d^/dx 5 ) 

= = = = .= .= _.=,.,_______ ) current scaled 5th derivatives 

D+6n-l d - 2 =m (h 751) d 5 y /dx 5 ) 
n ~~ ' n 


D+6n ) 

- = -_--_ = = „_„«.„__ ) guard digits for y. 

D+7n-l ) X 


D+7n ) 

----..-..--..-....»..- ) "predicted" values of 2~ m f . 
D+8n-l ) 1 


D+8n ) 

=.____- = _«,.,_.._„_«. ) temporary storage 

D+9n-l ) 


D+9n ) 

_--___-__._-__-■ ) values of y. at beginning of last 

D+10n-l ) step executed 
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The second and higher derivatives of y. alluded to in 
this description are to be interpreted as higher derivatives 
of an approximating polynomial for y. (x) and are developed 
automatically by the computer • The h which occurs in the 
description is the signed elementary increment to x which 
was used in the last elementary step done by the computer* 



AUXILIARY SUBROUTINES 

The auxiliary subroutine reads the values of x and y. from 

locations D-l to D+n-1 and computes the scaled derivatives 

—m 
2 f . and places them in locations D+n to EH-2n»i This 

subroutine must be supplied by the user since it depends 
on the particular system of differential equations he 
wants solved. It must be designed, for standard subroutine 
entry , with link address in the Q register, and to be 
entered by transferring control to the left side of loca- 
tion *5< . It must be capable of calculating the derivatives 
anywhere in the x-interval over which the equations are to 
be integrated, since the points at which the F7 routine 
will call for the derivatives are not known in advance. 
The accuracy of the auxiliary subroutine must be consistent 
with the accuracy demanded in the solution (a point some- 
times overlooked), and, this implies that 2~ f . must be 
calculated, correctly to e + m binary digits or 39 binary 
digits, whichever is the lesser number. 
INFORMATION AVAILABLE UPON EXIT FROM THIS ROUTINES 

Most of the useful information is contained in the data 
bank, and for such the reader is referred to the above 
description of the data bank. 

The "guard digit" locations D+6n+i-l are initially 
(automatically) set to l/2 so that the registers D+i-1 

will contain the rounded current values of the y. » If 

— — — — "'i 

the double precision value of y. is desired, say to improve 
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the accuracy of the auxiliary subroutine arithmetic, it 
may thus be found from the formula 

y. = N(D+i-l) + 2" 59 [N(D+6n+i-l) - l/2] 

The signed increment to x (signed elementary interval h) 
last used by the routine may be found in location (F7) 
+ 209, which is the last word of ¥J . . It is this value 
of h, including sign, which obtains in the data bank 
description. 

The recommended procedure for finding y. at an intermediate 
point j, i.e., a point interior to the elementary interval, 
is to use the available approximating polynomial to inter™ 

pOx3oe» mis is a c c omp-L 1 s x j.C'^i. as ioj l.owSo 

1.) Integrate until the point of interest lies beyond 

x-h but not beyond x. 
2 . ) Read out h and data bank locations D-l to D+6n-l 

inclusive, 
3«) Let the point of interest be x-9h| < 9 < 1. Then 

y. (x-9h) is given with accuracy equivalent to the 

accuracy of y. (x) , rf by the interpolation formula? 

y, (x-9h) = y, (x) - 2 m 9h [2~ m f , - 9 a, + 9 2 b. - oV+gV] 
1 1 11111 

Roots and other similar functions of the solution are also 

best found, by using the above interpolation procedure. 

Other procedures for treating intermediate points, such as 
continually integrating with smaller and smaller interval 
in such a way as to converge on the point of interest, are 
awkward, time-consuming, and inaccurate compared to the 
above recommended procedure. 



NOTES ON SCALING AND CHOICE OF m, JL , es 

o 



The parameter m is provided so that the y. and the f . 
may be scaled, independently for efficient use of the EH-i-1 and 
D-f-n+i-1 registers. An obvious requirement is that the 2 f., as 
well as the y. and x, must all be within register capacity 
throughout the range of integration. 
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In general, if less than full register accuracy is re- 
quired it is best to work near the less significant 
end of the registers , i.e., to scale the y. so that they 
are small. The advantage of this is that m may be taken 
small, which in turn permits Ji to be taken small if 
desired for increased speed. Round-off errors will not 
propagate more than a few places into the y. registers 
because of the guard digits kept. 

9 ■■ ''.should be chosen (a) large enough to guarantee 
-^ o 

availability of the solution of the printout or display 
intervals desired, and (b) so that 2 o is just larger 
than the largest elementary interval the computer will 
use, if the latter interval can be estimated,, The penalty 
for specifying JL too small is excessive automatic 
starting time, which is minor 3 whereas the penalty for 
J? too large is excess time spent by the computer in 
developing a more accurate solution than asked for. 

If some 2 f_ has a relatively sharp temporary fluctuation 

i 
in value in a narrow range of x (e.g., an approximation 

to a S-f unction) then one must guarantee, by taking JC 

large enough or otherwise, that at least one point 

x = x + (integer) 2" id falls within the region of the 

fluctuation! otherwise, there is some risk of the 

computer missing the fluctuation entirely.. 

Discontinuities or "jumps" in 2 f ., if any, must not 
exceed l/8 in magnitude so that the higher-derivative 
registers will not overflow during the process of Integra- 
tin past the discontinuity . 

e should be chosen after the scaling for the y. has been 
chosen. If the maximum value of y. is 2 then there 
will be at least about (e-*-d) correct significant bits 
of the y. developed after integrating from x to x . + 1. 
The maximum error tends generally to increase linearly with 
the range of integration, except that the number of correct 
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significant digits cannot be expected to exceed the 
initial number of significant digits. For 'Veil -behaved" 
systems of equations s aamely ones such that 



-m 



df, 



dx 



= ^1 



the error is generally less than 2 s typically 2 
For very ill -behaved equations ^ such that 



-m 



df, 
dx 



^2 56 



the error will in general be larger 3 typically 2 



■e+1 



A more precise estimate of the error in a particular 
result may be got by using two different e's in turn and 
observing which digits in y. change • e should be changed 
by at least 3 or k units in this process « 

The relationship between the value of e specified, and the 
interval automatically chosen by the computer is as 
follows. The interval h chosen is the largest interval 
satisfying all of the three conditions % 



h < 2 



-K. 



h 



1>f. 



< 0.38| h |d y./dx I < 3 • 2' 



In the last of these conditions the 6th derivative is to 
be interpreted as the 6th derivative of an approximating 
polynomial developed by the computer for y. (x). The 
second of the above conditions guarantees the stability of 
the numerical procedure and minimizes inherited errors 
in cases where such errors are strongly magnified . The 
third condition, :bounds the truncation error in such a 
way as to supply the final accuracy specified via the 
accuracy parameter e. 



- 9 - 

TO FREEZE THE INTERVAL; 

The following modifications will transform the F7 routine 
into a fixed-interval routine operating with elementary 
interval h = + h = + 2 ^o. The accuracy parameter e 
and all interval controlling tests will be ignored. 
At; Overwrite: With? 

(F7) + kQ 32 55L F5 139L % ((F7) + 7*0 F F5 ((F7) +■ 139) F 

(F7) + 81 ^0 110L hO 207L 50 F ^0 ((F7) + 207) F 

Warning : The accuracy and stability of the procedure 
may suffer in such a fixed- interval application if the 
condition 



h ~i 
°^3 



< o38 



is not satisfied. 



DESCRIPTION OF ROUTINE: 

A detailed description of the method used in this routine 
may he found in University of Illinois Coordinated Science 
Laboratory Report R-127, "On Numerical Integration of 
Ordinary Differential Equations, "'May, 19&L, "by Arnold 
Nordsieckj also to be published, in "Mathematics of 
Computation," Journal of the Division of Mathematics, 
National Academy of Sciences - National Research Council, 
Washington, D. Co, about October, 1961* 

The method works with polynomials of 5th degree approxima- 
ting to each of the y., and these polynomials are specified 
by their Oth to 5"th derivatives in order to facilitate 

changes of interval. (Actually the derivatives beyond the 

— m 
first are multiplied by 2 and by a power of the current 

elementary interval h and divided by a factorial, as 

indicated in the data bank description, in order to keep 

them in machine range and to facilitate the arithmetic). 

Equations with discontinuous derivatives are integrable 
by this method with accuracy comparable with the accuracy 
for the case of continuous derivatives provided only that 
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the jumps in the derivatives are finite and are scaled 
according to the restriction stated on page 7° 

The elementary step from>x to (x + h) consists of applying 
the working equations: 

y (x + h) = y.. (x) . 

1 -L - 

+ 2 m h (2" m f.(x) + a.(x) + b.(x) + c.(x) + d.(x) + 2^ [ 2 " m f. (x + h)=-2" m f . P ] ) 



2° m f .' = 2" m f . (x) + 2a. (x) + 3b. (x) + 4c„ (x) + 5d. (x) 



2p r ~m^ / ,, x ~m^ Pi 



a.(x + h) = a.(x) + 3b. (x) + 6c. (x) + 10d.(x) + f§[2~ m f. (x+h)-2~ m f .*] 
b^x + h) = b i (x) + ^ ± (x) + 10d i (x) + ||[ " ] 

c. (x:-+ h) = c.(x) + 5d ± (x) + ^[ " ] 

d.(x + h)= \w+h [ " j 

Here the quantities evaluated at x are known from the 
: previous elementary step (except in the -very first step,' 
in which case the "initial entry" must he used and an 
automatic starting procedure is invoked! see "below )$ and 
the quantities at (x + h) are to be found. The coefficients 
95/288, 23/2^ 35/72, 5/1*8, 1/120 are specially chosen 
to maximize the stability of the method, to optimize the 
process of integration across finite discontinuities of 
the derivatives, and to make the truncation error O(h') 
rather than 0(h ) as it would be for arbitrary coefficients . 
The truncation error in y. per elementary step is 

dx 
provided the seventh derivative exists 

The solution of the working equations proceeds in three 
stages. Stage 1 consists of "predicting" all six quantities 
y. ... d. at x + h, i e , applying the working equations 
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without the [ ] ■. terms , using as tentative value of h the 
value which was accepted, in the last previous step, or 
twice that value if the conditions for doubling h were 
fulfilled. Stage 2 consists of solving the complete first 
working equation, which is an implicit equation since the 
f . (x + h) depend .on the y. (x + h) in virtue of the differen- 
tial equations o The equation is Solved iteratively by 
inserting the "predicted" y. on the right, appealing to 
the auxiliary subroutine for f . (x + h), thus producing 
first improved y. (x + h), appealing again to the auxiliary 
subroutine and thus finally getting second improved 
y„ (x + h). The iterative procedure is always terminated, 
after just these two iterations, the second improved, 
y. (x + h) and. f .(x + h) being accepted, as adequate approxi- 
mations . 

At completion of stage 2 two tests are made to determine 
whether h is sufficiently small. One test determines 
whether or not the iterative solution of the implicit 
equation was sufficiently convergent, and. this comes to 
whether or not 

Tff. 



95 
2BB 7 






< 1/8 



because the left side of this inequality is the conver- 
gence factor* The second test determines whether the 
truncation error is small enough to be consistent with 
the accuracy parameter e, and this comes to whether 

|f ± (x + h) - f ± p | < (3^03) 2~7Tn| 

If either of these tests is violated we discard the 
computations of stages 1 and 2, halve h and enter upon 
stage 1 again. If both tests are satisfied we ascertain 
whether they are 'Wersatisf ied, " i.e., whether a doubled 
h would likely satisfy them, and if so note this fact 
for future reference; then proceed to stage 3> which 
consists of "correcting" a., b„, c, d. by adding the 
[ ] terms » 
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The above discussion indicates the conditions for 
automatic reduction of interval* If the tests are not 
satisfied even after the interval has been reduced to 
2 the FF stop is invoked. The conditions for increasing 

the interval, on the other hand, are as follows; If 

i . -1 
before entering stage 1 we have |hj < 2 o and both 

tests were previously "oversatisfied" and the digits 

of (x-x ) are such that doubling h will not cause the 

-Q 

next point x = x + (integer) 2 ^"o to be missed,, then 
h is doubled before stage 1 is entered; but after a 
reversal of sense of integration h may not be doubled 
until k elementary steps have elapsed. The last condition 
is necessary to discourage occasional erratic interval 
behavior after reversal e 

Automatic starting is a special procedure automatically 
invoked before the first step away from the initial value 
x , when neither the correct h nor the quantities a.. . . • d. 

o 11 

are known* This feature relieves the user of the task of 

supplying special starting information and. requires 

him to supply only the logically essential initial values 

of the y. . 
J i 

In the automatic starting mode the computer clears the 

a. ... d. locations and then repeatedly integrates four 

steps forward, and. four steps back to x « Because of 

the high degree of stability of the method the quantities 

a. . •• d. converge rapidly to the correct values in this 

process 5 which is thus a successive approximation method. 

of a special sort for fitting a 5th degree polynomial 

to y. (x) near x . During starting the interval is also 

decreased as necessary, so that when the process is 

complete the correct h as well as the correct a. °oo d. 

have been established. Thereafter one meaningful forward. 

step is taken and control is returned to the link address « 

From the point of view of the user the first step is 

just like any other step except that a) he must supply 
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the initial conditions and the program parameter and 
enter by the "initial entry, " and b ) it takes 2h or 
more times as long as other steps. 

Complete details of the starting process and of the whole 
method in general are given in the paper "On Numerical 
Integration of Ordinary Differential Equations" referred 
to above. 
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00K(F7) 
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00 F 




Interlude 




00 L 
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80 IF 




Construct constants 
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FF 77F 
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1 


49 136L 
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K5 F 


by 2 


£— reversing entry 


2 


l*-2 1L 

Fl 207L 




/ -.save link address 


3 


32 i6l 

L5 F 




transfer unless initial entry 






k 
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reset JL- 




5 


10 12F 








6 
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reset X~- m + 1 to 

* program 
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7 


00 F 
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9 


09 IF 










10 F 


by 6' 


address = e + 1 

address = X* - m + 1 





10 


00 F 


by 6 




11 


4o 208L 
L5 1023S 3 




reset test quantity t ________ 


1 




40 1022S3 




I save x 
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17 
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undo prediction 


62 


50 91L 






F5 62L 






63 


26 27L 
L5 208L 




J 




/ 


S 


6^4- 


80 IF 
32 65L 




2t -> t with ceiling 1 




65 


40 208L 
F5 209L 




\ 




66 


10 IF 
h-0 209L 




> l/2 Ax ->Ax and FF stop 




67 


00 40F 




if underflow halve h 




36 L 


. - 


/ 




68 


L5 7L 
Lk 199L 


' 


\ JL.^Jl 




69 


^0 7L 


J 






70 


L5 198L 
22 26L 
L5 20 7L 




l/2a ± ^a ± , l/bb ->b. etc. 




71 


32 72L 
LI 209L 




transfer on l6th starting step 












— —— 



LOCATION 



72 

73 
7^ 

75 
76 
77 
78 
79 
80 
81 
82 

83 

8^ 

85 
86 

87 
88 

89 



ORDER 



22 39L 
L5 H3L 
32 13L 
22 15L 
50 118L 
F5 7^L 
26 27L 
L5 207L 
80 IF 
kO 136L 
kO 207L 
36 39L 
80 IF 
kO 207L 
32 17L 
80 IF 
32 81L 
80 IF 
kO 110L 
^0 207L 
50 86L 
F5 82L 
26 27L 
L5 110L 
32 63L 
22 15L 
IA 205L 
26 87L 
LO 205L 
k6 90L 
ko 88l 

11 IF 
L4 F 
kO F 
L5 207L 
32 30L 



NOTES 
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F 7 



adjust x 



transfer if truncation error too large 

"correct" a., b., c., d. 

beginning of starting mode control 
to indicate "not reverse" 
transfer if code digit is 

transfer if 10 

transfer if 110 

to distinguish between step l6 and step 2k 



reinsert correct initial y; 



transfer if step 16 

end of starting mode control "—"" ^ 

^ entry for "save" 



} 



by 22', 38', 
87 



£ — entry for "reinsert" 



clear accumulator 



"save y °* 
^i 

or "reinsert 

process 
subroutine 



transfer if starting mode 



LOCATION 


ORDER J 




NOTES PAGE 7 F? 


90 


k9 F ' 


by 86 ' 


set guard, digits = 1/2 




91 


22 30L 
^0 111L 




\/ 




/N 






LO 20^L 








92 


k6 105L 
LO 203L 








93 


k6 lO^L 
k2 lk6L 








9k 


L^ 20^L 
k6 107L 








95 


^2 107L 










40 112L 




"prediction" 


96 


L^ 20*kL 




process 




k6 108L 




subroutine 


97 


k-2 108L 




(through 




lj-0 110L 




117) 




98 


L4 S03L 
\2 109L 








99 


L^ 20^L 
k2 105L 








100 


^0 113L 
L^ 203L 








101 


\2 ll^L 
L4 203L 








102 


4-2 104L 

ii-6 i46l 








103 


L5 133L 
36 107L 




transfer if undoing prediction 




10^ 


L5 F 


by 93 


} 






ho F 


by 102 


1 save y ± 




105 


L5 F 
40 F 


by 92 
by 99' 


and f . 

J 




106 


L5 136L 
32 109L 




transfer unless reversing 




107 


LI F 


by 9^' 

j 


1 






kO F 


"by 95 


)r -a . -> a . and -c . -> c . 

11 11 





' LOCATION 


ORDER 


1 


NOTES PAGE 8 


F 7 


108 


LI F 
40 F 


by 96' 
lay 97 


J 






109 


10 IF 




waste order 








L5 F 


by 98' 


~ 






110 


L4 F 
40 F 


| by 81 




"prediction 


111 


L4 F 


1 by 60,91 




■^process 




40 F 


J 50,177' 


predict first y., then 


subroutine 


112 


L4 F 


1 by 49S95S 


f . , ... then c . 








4o F 


J 1^9 


r 






113 


L4 F 
40 F 


| by 100,147' , 
J 175 








ll4 


22 ll4L 


by 117 


alternate address l4^L 






115 


^4-0 F 

10 F 
00 IF 


by 101 
by 14' 


address =Jc - m + 1 






116 


4o IF 
L5 20L 




; 






117 


46 ll4L 
22 109L 






\jz 


T 


118 


4o 136L 












LO 203L 






"correction" 


119 


1+0 139L 

L4 204L 




■ 


process 
subroutine 


120 


40 133L 
L4 204L 










121 


te l4lL 
LO 20 3L 










122 


4o 131L 
LO 205L 






"correction" 
process 


123 


46 l4lL 
LO 203L 






subroutine 
(through 


124 


46 125L 






145 left) 




42 125L 










125 


LI F . 


by 124 










L4 F 


by 124' 






1 "*" 



LOCATION 
126 

127 
128 

129 
130 
131 
132 
133 

135 
136 

137 
138 

139 

1^1 
1^2 
1^3 



ORDER 



NOTES 



PAGE 9 

■P:. 



F 7 



"H 



by 122 



by 61' , 
120 



2" m [f.(x+h) - f. 1 : -> location 



transfer if correction terms 5 



transfer if correction terms > 
l/2 - l6/l20, approxo, -> accumulator 

[2" m (f-fP); + !6 x 2~ 39 ] + 1/120 for d ± 
for later 75 order 



d. corrected 

1 



c . corrected 
1 



waste order; left function used 



L "by 1, 76', 
118 



b . corrected 

1 



22 30L 

F5 IF 



by k5 % ,119, 
183 



by 123 
by 121 



a . corrected 

1 



transfer if in starting mode 



adjust guard, digits "correction" 

process 
insert exit link address subroutine 



LOCATION , 


ORDER 




NOTES PAGE 10 


F 7 


144 


kO IF 
LJ 134L 




l/2 + 16/120 approx, 
-> accumulator 






lij-5 


22 129L 




t ~\ "l *n • _o • i _____ _i_ * 


\J/ 




L5 IF 


"by 20 


LI IF if integrating 


backwards 

•'K 


146 


Lk F 


by 102' 










40 F 


by 93 f 


y. predicted 






147 


L5 151L 






part of 




40 113L 




blocking order -» 113 


"prediction"! 


ikQ 


L5 117L 
26 117L 






process 
subroutine 


149 


40 112L 
22 109L 




blocking order -» 112 






150 


kO 111L 
22 109L 




blocking order -» 111 






1 Rl 


L5 152L 




^ 






J -^'- J - 












26 l49L 




* constants for 






152 


L5 1^3L 




blocking 








26 15.OL 




/ 


\U 


/K 


153 


LO 205L 
k6 166L 










15^ 


42 I78L 
42 179L 










155 


LO 203L 












46 160L 






"iteration" 


156 


L4 204L 
46 161L 






process 
subroutine 


157 


42 l6lL 
L4 203L 






(through 
183) 


158 


46 165L 
46 171L 










159 


L4 203L 
46 178L 










160 


50 F 

00 8f 


by 155' 








l6l 


LI F 


by 156' 










L4 f 


isy 157 









LOCATION 


ORDER j 
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162 


kO F 

Sk F | 




2" m [f.(x+h) - f. p ] ^location 
i x 1 




163 


kO 188LJ 
50 F 




for oscillation test 




161* 


7J 191L 
kO F 








165 


L^ F 
kO IF 


by 158 






166 


51 F 


"by 153' 








10 IF 




"iteration" 


167 


00 F 


by k2 


address = JC ~ m + 1 process 




■lA 1L 




LO IF if going backwards subroutine 


168 


10 F 
00 IF 


by k2 l 


address = \_ - m + 1 




169 


ko IF 
Fl 112L 




Ay* —> location 1 




170 


36 I78L 

S5 F 




transfer if 1st iteration 




171 


kO F 
L7 F 


by 158' 


save new guard digits 




172 


80 IF 
L2 188L 








173 


36 <?kL 

L3 F 




transfer if oscillation 




17 1 * 


LO 208L 
36 I76L 


i 


fsmall enough 
transfer if truncation error 




175 


kl 113L 
26 I78L 




note test failure 




176 


10 6f 

L2 F 








177 


36 I78L 




transfer if test oversatisfied. § 




kl 111L 




veto h-doubling 




178 


L5 F 


by 159' 








50 F 


by 15k 




, 


179 


L^ IF 










kO F 


by 15 V 


y. (x+h), 1st or 2nd iterate f 


.... 

u__ 




1 
■4 


1 im fr 



LOCATION 



180 
181 
182 
183 
iQk 
185 
186 
187 
188 
189 
190 
191 
192 
193 
19^ 
195 
196 
197 



ORDER 



NOTES 
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max change 
in y\ -» 139L 



"interation" 

process 

subroutine 



\V 



'TV 



"modify 



3 <V 



d. " process 
subroutine 



by 186 ' 

by 26' ,163, alternate functions 10, 00 

18V,185,18?'/ 



transfer unless d. done 



V 



00 io*fi 6666 6667 J 

00 f I 

00 ^861 1111 liiu 

00 F I 

00 ki6 6666 666 7J 

11 l8^L 

ko 15L 

00 l8^L 

ko 50L 



for clearing a. 



for doubling h 



d. 

1 




NOTES 

ii- for halting h 
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F 7 



increment 



for changing forward ^— •£ backward. 



starting code word 

causes SADOI to put following tape 
material immediately after this program 



end of tape 



constants depending 
> on n and D, constructed 

by interlude 



starting code and 
start/run switch 
-e + *~ -m 



-t - -2 



,JL 



Ax z. h = + 2 



end of F7 program 



